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achieved  to  date,  especially  for  the  case  of  spin  models  with  short-range  couplings.  However,  progress 
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this  work  must  maintain  toward  the  development  of  a  comparable  understanding  in  long-range  interacting  models,  m 
author (1°  and  the  title  of  particular  out-of- equilibrium,  remains  limited.  In  a  recent  work,  we  proposed  a  semiclassical 
the^woirk,  journal eitation  numerical  method  to  study  spin  models,  the  discrete  truncated  Wigner  approximation  (DTWA),  and 
demonstrated  its  capability  to  correctly  capture  the  dynamics  of  one-  and  two-point  correlations  in 
one-dimensional  (ID)  systems.  Here  we  go  one  step  forward  and  use  the  DTWA  method  to  study  the 
dynamics  of  correlations  in  two-dimensional  (2D)  systems  with  many  spins  and  different  types  of 
long-range  couplings,  in  regimes  where  other  numerical  methods  are  generally  unreliable.  We 
compute  spatial  and  time-dependent  correlations  for  spin- couplings  that  decay  with  distance  as  a 
power-law  and  determine  the  velocity  at  which  correlations  propagate  through  the  system.  Sharp 
changes  in  the  behavior  of  those  velocities  are  found  as  a  function  of  the  power-law  decay  exponent. 
Our  predictions  are  relevant  for  a  broad  range  of  systems  including  solid  state  materials,  atom-photon 
systems  and  ultracold  gases  of  polar  molecules,  trapped  ions,  Rydberg,  and  magnetic  atoms.  We 
validate  the  DTWA  predictions  for  small  2D  systems  and  ID  systems,  but  ultimately,  in  the  spirt  of 
quantum  simulation,  experiments  will  be  needed  to  confirm  our  predictions  for  large  2D  systems. 


1.  Introduction 

An  important  advance  toward  understanding  non -equilibrium  phenomena  has  been  made  possible  by  recent 
advances  in  cooling,  trapping  and  manipulating  atomic,  molecular  and  optical  (AMO)  systems  [  1  ] .  In  contrast 
to  solid  state  systems,  where  studying  equilibrium  situations  is  the  default  approach  (given  their  complex 
environment  and  fast  relaxation  rates),  AMO  systems  provide  a  unique  platform  to  observe  and  investigate  non¬ 
equilibrium  quantum  dynamics  in  strongly  interacting  many-body  systems  [2] .  Their  high  dynamic  tunability 
even  during  the  course  of  an  experiment,  their  decoupling  from  the  external  environment  and  their 
characteristic  low-energy  scales,  lead  to  long  non- equilibrium  time- scales  over  which  the  system  can  be  followed 
almost  in  real  time.  Quantum  quenches,  i.e.  the  dynamics  induced  by  abruptly  changing  parameters  of  the 
system,  is  currently  a  common  protocol  used  to  probe  AMO  systems. 

One  of  the  most  promising  opportunities  offered  by  modern  AMO  physics  is  the  ability  to  engineer 
interatomic  interactions  different  from  the  standard  contact  and  isotropic  interactions  arising  from  ultracold 
collisions.  At  the  heart  of  this  capability  are  recent  experimental  developments  on  controlling  AMO  systems 
with  complex  internal  structure  and  with  enlarged  sets  of  degrees  of  freedom  such  as  polar  molecules  [3] , 
trapped  ions  [4-6] ,  magnetic  atoms  [7-9] ,  Rydberg  atoms  [  10-15] ,  and  alkaline  Earth  atoms  [16-20] .  All  these 
systems  have  in  common  that  they  can  exhibit  long-range  interactions.  This  experimental  progress  is  opening 
new  frontiers,  and  at  the  same  time  demanding  for  improved  theoretical  techniques,  that  are  capable  of  dealing 
with  the  complicated  non-equilibrium  quantum  dynamics  of  long-range  interacting  systems. 
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In  a  prior  work  [2 1  ] ,  we  proposed  a  semiclassical  phase-space  method  to  study  non -equilibrium  quantum 
dynamics.  We  used  this  numerical  method,  that  we  named  the  discrete  truncated  Wigner  approximation 
(DTWA),  to  study  the  dynamics  of  single  particle  observables  and  correlation  functions  after  a  quench.  The 
DTWA  was  benchmarked  in  one- dimensional  (ID)  spin  models  with  numerically  exact  time- dependent  density 
matrix  renormalization  group  (t-DMRG)  calculations  [22-25]  and  excellent  agreement  was  found. 

In  this  work,  we  generalize  the  calculations  to  two-dimensional  (2D)  Ising  and  XT  spin  models  with  various 
ranges  of  interactions.  We  study  the  time- evolution  in  a  setup  that  is  equivalent  to  a  Ramsey- type  procedure,  as 
realized  in  recent  experiments.  This  dynamical  protocol  has  been  used,  for  example,  to  observe  dipolar  spin- 
exchange  interactions  in  ultracold  molecules  [3, 26, 27] ,  to  benchmark  the  Ising  dynamics  of  hundreds  of 
trapped  ions  [5] ,  and  to  precisely  measure  atomic  transitions  as  well  as  many-body  interactions  in  optical  lattice 
clocks  [  1 7, 1 8, 28, 29] .  Using  the  DTWA  we  compute  the  dynamics  of  the  collective  spin  as  well  as  spatially 
resolved  two-point  correlation  functions.  Since  the  applicability  of  the  t-DMRG  method  becomes  limited  in  2D, 
to  benchmark  the  DTWA  we  perform  numerical  comparisons  with  small  systems  (where  exact  diagonalization 
(ED)  is  possible),  and  with  the  analytically  solvable  Ising  case  [30-32] .  We  then  extend  the  calculations  to  large 
XT  spin  models,  and  find  remarkably  sharp  changes  in  the  propagation  of  correlations  as  we  vary  the  power  law- 
decay  exponent  of  the  interactions. 

The  dynamics  of  the  two-point  correlations  is  directly  linked  to  the  speed  of  propagation  of  information  in 
quantum  many-body  systems,  a  topic  of  great  interest  to  quantum  information  science  and  currently  subjected 
to  intensive  investigation.  For  systems  with  short-range  interactions,  there  is  a  well  understood  bound  (derived 
by  Lieb  and  Robinson)  that  limits  correlations  to  remain  within  a  linear  effective  ‘light  cone’  region  [33, 34] .  On 
the  contrary  there  are  many  open  questions  about  what  limits  the  propagation  of  information  in  quantum 
many-body  systems  with  long-range  interactions  [3 5-45] .  Here  we  use  the  DTWA  method  to  determine  the 
shape  of  the  causal  region  and  the  speed  at  which  correlations  propagate  after  a  global  quench.  We  compute  the 
full  crossover  of  the  dynamics  when  changing  the  range  of  the  interactions  over  a  large  range  in  ID  and  2D 
systems,  and  observe  remarkable  agreement  of  DTWA  results  with  t-DMRG  predictions  in  the  ID  case.  Our 
calculations  and  their  natural  variations  (e.g.  local  instead  of  global  quenches)  should  be  testable  in  experiments 
with  polar  molecules  and  trapped  ions  in  the  near  future. 

This  paper  is  organized  as  follows:  in  section  2  we  introduce  the  models  that  are  studied.  In  section  3  we 
review  the  DTWA  technique  that  was  introduced  in  [2 1  ].  We  benchmark  the  DTWA  method  by  comparing  the 
dynamics  of  single-spin  observables  and  correlation  functions  with  exact  solutions  in  section  4.  In  section  5 
those  calculations  are  extended  to  large  systems  where  currently  no  other  method  is  applicable.  In  section  6  we 
use  the  DTWA  for  a  systematic  calculation  of  the  light- cone  dynamics  as  we  vary  the  range  of  the  interactions. 
Finally,  section  7  concludes  and  provides  an  outlook. 

2.  Spin  models  and  dynamics 


We  will  focus  our  attention  on  Hamiltonians  that  fall  under  the  generic  heading  of  spin- 1/2  XXZ  models  given 
by  [fi  —  1) 


(1) 


where  the  sum  extends  over  all  pairs  of  sites  of  an  arbitrary  lattice,  6?'y'z  are  Pauli  matrices  for  the  spin  on  site  z, 

Jij  =  Jji& nd  Ja  =  0.  In  our  analysis  the  interactions  are  assumed  to  decay  as  a  function  of  the  distance  with  a 
decay  exponent  a,  that  is  Jtf’z  =  J  (a  I  |  |  )a.  Here,  is  the  vector  connecting  spins  on  sites  z  and  j  and  a  is  the 

lattice  spacing.  We  concentrate  our  study  on  two  specific  cases,  Ising  (Jtf  =  0)andXT  (/|  =  0)  interactions.  We 
consider  a  general  2D  grid,  i.e.  a  lattice  with  Nx  X  Ny  =  M  sites  with  spins  at  positions  r \  =  a  (nXi  ny ),  where 
nXyy  are  integers.  The  lattice  spacing  a  is  set  to  1  throughout  this  paper. 

Spin- 1/2  XXZ  models  broadly  describe  a  variety  of  physical  systems.  For  instance,  in  the  AMO  context,  XXZ 
spin  Hamiltonians  have  been  used  to  model  the  dynamics  of  ultracold  molecules  in  optical  lattices,  trapped  ions, 
Rydberg  atoms,  neutral  atoms  in  optical  clocks,  and  ultracold  magnetic  atoms.  A  summary  of  how  these  models 
are  realized  in  those  systems  can  be  found  in  [46] .  Here  we  describe  the  two  most  relevant  ones  for  this  work: 
ultracold  polar  molecules  and  trapped  ions. 

In  ultracold  polar  molecules  pinned  in  optical  lattices,  the  spin- 1/2  degree  of  freedom  can  be  encoded  in  two 
rotational  states,  and  the  spin-spin  couplings  are  generated  by  dipolar  interactions.  The  difference  in  dipole 
moments  between  the  two  states  (which  arises  in  the  presence  of  an  electric  field)  generates  the  Ising  term  while 
transition  dipole  moments  between  the  two  rotational  states  (which  can  exist  even  in  the  absence  of  an  electric 
field)  give  rise  to  the  spin-exchange  terms  [47] .  The  ratio  between  the  Ising  and  XT  couplings  can  be 
manipulated  using  electromagnetic  fields  [48-53] .  The  case  without  electric  field  which  implements  the  pure  XT 
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model  has  been  recently  realized  with  KRb  polar  molecules  in  a  3D  optical  lattice  [3, 27] .  In  general,  dipolar 
interactions  are  long-ranged  and  spatially  anisotropic.  In  a  2D  geometry,  however,  they  become  isotropic  if  the 
electric  field  that  sets  the  quantization  axis  is  set  perpendicular  to  the  plane  containing  the  molecules.  This  is  the 
case  considered  throughout  this  paper. 

Crystals  of  2D  self- assembled  trapped  ions  can  also  be  used  to  implement  specific  cases  of  equation  ( 1 ),  see  [5] . 
By  addressing  the  ions  confined  by  a  Penning  trap  with  a  spin- dependent  optical  potential,  the  vibrations  of  the 
crystal  mediate  a  long-range  Ising  interaction  that  can  be  approximately  described  by  a  power-law  with 
0  ^  a  <  3  [4, 5, 54—57] .  To  engineer  an  XY model,  one  needs  to  add  a  strong  transverse  field  that  projects  out 
the  off- resonant  terms  in  the  Ising  interactions  that  change  the  magnetization  along  the  field  quantization 
direction  [40,41]. 

The  dynamical  procedure  considered  here,  is  identical  to  a  Ramsey  spectroscopy  setup.  It  has  been 
implemented  in  various  recent  experiments  as  a  diagnostic  tool  for  interactions  [27] .  It  consists  of  preparing  an 
initial  state  with  all  spins  aligned  (at  time  t  =  0)  along  a  specific  direction,  here  we  consider  it  to  be  the  x  direction, 
i.e.  \y/(t=  0))  =  0-^  (|  t  )i  +  I  i  )i)/V2.  Then  this  initial  state  evolves  under  the  Ising  or  XTHamiltonian 
( 1 )  for  a  time  t,  leading  to  the  state  |  \j/  (t) )  =  exp  (— itH)  \  y/  ( t  =  0) ).  Afterwards  one  measures  expectation 
values  of  an  observable  with  the  time-evolved  state,  (y/(t)  \  O  \y/  (t) ).  In  this  paper,  we  focus  on  two-point 
correlations  and  the  collective  spin  along  x  as  observables. 

3.  The  DTWA  method 

Phase  space  methods,  such  as  the  truncated  Wigner  approximation,  solve  the  quantum  dynamics  approximately 
by  replacing  the  time -evolution  by  a  semi- classical  evolution  via  classical  trajectories.  The  quantum  uncertainty 
in  the  initial  state  is  accounted  for  by  an  average  over  different  initial  conditions  [58, 59] ,  determined  by  the 
Wigner  function.  Although  the  Truncated  Wigner  Approximation  approximation  was  initially  developed  to 
deal  with  systems  with  continuous  degrees  of  freedom  [60, 61  ] ,  it  has  also  been  adopted  to  treat  collective  spin 
models.  In  this  case  the  standard  method  has  been  to  approximate  the  Wigner  function  by  a  continuous 
Gaussian  distribution  that  facilitates  the  sampling  of  trajectories  [59] .  This  continuous  approximation,  however, 
misses  important  aspects  inherent  to  the  discrete  nature  of  spin  variables  and  it  is  unsuitable  in  systems  with 
finite-range  interactions.  To  deal  with  more  generic  types  of  spin  models,  recently  we  proposed  instead  to 
sample  a  discrete  Wigner  function  for  each  spin  and  named  this  approach  the  DTWA  method  [2 1  ] .  In  this 
section  we  present  an  overview  of  the  DTWA  method.  For  details  the  reader  is  referred  to  [2 1  ] . 

Operators  and  wave  functions  on  the  Hilbert  space  of  a  quantum  system  can  be,  equivalently,  represented 
(mapped)  on  a  classical  phase  space.  There,  any  operator  O  corresponds  to  a  real-valued  function  of  the  classical 
phase-space  variables,  Ow  (a  so-called  Weyl  symbol).  The  phase-space  function  corresponding  to  the  density 
matrix  is  precisely  the  Wigner  function  w.  For  continuous  variables  p,  q  in  one  dimension  (for  simplicity  of 
presentation),  the  expectation  value  of  an  operator  can  be  exactly  represented  as 
<6)(f)  =  J fdpdq  w(p,  q;  t)Ow  (p,  q). 

For  quantum  systems  with  discrete  degrees  of  freedom,  one  can  introduce  a  'discrete  phase- space’  in  various 
ways,  see  [62]  and  references  therein.  Here  we  use  the  representation  of  Wootters  [62, 63] .  For  a  single  spin- 1/2, 
it  uses  four  distinct  phase-points,  and  all  phase-space  functions  are  thus  defined  as  2  X  2  matrices.  In  our 
approximation,  the  phase-space  of  iV spins  factorizes  into  a  product  of  N phase- spaces  for  each  individual  spin. 

In  both  continuous  and  discrete  variables  however,  it  is  not  possible  to  compute  the  time- evolution  exactly. 
The  spirit  of  the  truncated  Wigner  approximation  and  the  DTWA  is  to  take  quantum  fluctuations  into  account 
only  to  lowest  order  [59] .  In  particular,  we  switch  to  a  'Heisenberg  picture’,  such  that  the  Wigner  function  does 
not  evolve  in  time  (i.e.  it  is  fixed  to  the  initial  state)  while  the  operator- functions  are  time-dependent.  The 
approximation  that  we  make  is  to  assume  that  the  operators  in  phase  space  follow  their  classical  evolution: 

<o>(f)  =  2>(r;  o )ow(r;  t)  *  2 w(r,  o )ow’d-(r;  t),  (2) 

7  7 

where  y  runs  over  the  points  of  the  discrete  phase  space,  w  (y;  0)  is  the  Wigner  function  at  t  =  0  on  the  discrete 
many-body  phase  space,  and  Ow,cl  (y;  t)  is  the  classically  evolved  operator- function  (Weyl  symbol)  that 
corresponds  to  our  observable. 

Equation  (2)  is  solved  numerically  by  choosing  a  large  number  nt  of  random  initial  spin-configurations,  with 
probability  according  to  w  (y;  0).  Each  of  this  'Monte-Carlo  trajectories’  is  evolved  independently  following  the 
classical  equations  of  motion  (see  below).  The  expectation  value  in  equation  (2)  is  calculated  by  averaging.  We 
find  that  the  number  of  required  trajectories,  nt ,  does  not  depend  on  the  system  size,  but  rather  on  the  observable 
under  consideration. 
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To  apply  the  truncated  Wigner  approximation  we  have  to  compute  the  classical  equations  of  motion  for  the 
spin  components  of  each  spin  i :  sf,  sf>  sf .  One  way  to  do  this  is  to  replace  spin  operators  3  by  classical  variables 

(s  =  (3))  in  the  Hamiltonian  and  to  compute  the  Poisson  bracket  [59] ,  giving  sf  =  2  Y  espysf  with  e  the 

dsf 

fully  antisymmetric  tensor.  Alternatively,  the  same  classical  (mean- field)  equations  of  motion  can  be  obtained 
from  a  product  state  ansatz  of  the  density  matrix  [21].  For  the  Ising  interaction  Hamiltonian  the  classical 
equations  for  the  spin  components  are  given  by: 


Sx  =  -2  sy  yp  sz  =  -2sy6z 

m 

=  2  ^Jn,msm  =  2 


(3) 

(4) 


SZn  =  0,  (5) 

where  we  introduced  the  quantity  /?'5=x’y,z  =  y  /«,m5m=x’r’z  which  can  be  interpreted  as  an  effective  magnetic 
mean-field  acting  on  spin  n  induced  by  the  other  spins.  For  the  XT  interaction  the  classical  equations  of  motion 
are  given  by: 


sx  =  2sz  Yf1  sr  =  2  sz6y 

m 

sy  =  ~2sz  Yl1  =  -2  szfix 

Jn  ^  Jn  ^  n  rn  > 

m 

sz  =  l'VlL  (sxsy  -  sysx)  =  2syBx  -  sx6y 

n  ^  \  JmJn  JmJn  )  —  ^  Jn  rn  Jn  rn  • 


(6) 

(7) 

(8) 


Note  that  the  sums  exclude  the  term  m  =  n  since  we  set  Jnn  =  0. 

In  practice,  applying  the  DTWA  means  to  solve  these  equations  of  motion  nt  times,  while  each  time  choosing 
a  different  random  initial  configuration.  For  our  particular  initial  state  (pointing  along  the  x  direction),  the 
prescription  for  the  correct  sampling  is  to  randomly  pick  values  of  the  orthogonal  spin-components  for  each 
spin  from  s^(0),  sz( 0)  E  {  -1,  +1}. 

The  error  of  the  DTWA  method,  i.e.  the  deviation  from  the  exact  solution  arises  entirely  due  to  the 
semiclassical  approximation  for  the  time -evolution  (cf  equation  (2)).  The  Wigner  function  of  the  initial  state,  on 
the  other  hand,  is  sampled  exactly  here  up  to  statistical  errors  (which  can  be  controlled  by  increrasing  the 
number  of  trajectories).  For  the  exactly  solvable  Ising  model,  the  DTWA  method  turns  out  to  be  able  to 
reproduce  the  exact  solution  for  single- particle  observables;  for  two-particle  correlations  the  error  can  be  given 
explicitly  [21]  (see  below) . 

We  note  that  the  DTWA  clearly  goes  beyond  the  mean-field  predictions.  A  pure  mean-field  theory  is  not 
only  incapable  to  capture  spin-spin  correlations  (they  are  all  zero  due  to  the  factorization  approximation  of  the 
density  matrix)  but  even  the  single  particle  mean-field  obsevables  can  be  completely  incorrect.  For  example 
equations  (6)-(8)  predict  no  dynamics  at  all  in  our  Ramsey  setup  where  the  collective  Bloch  vector  points 
initially  along  x. 


4.  Benchmarking  the  DTWA 

4.1.  Contrast 

Before  discussing  results  obtained  by  the  DTWA  method  for  the  propagation  of  correlations  through  a  2D 
system,  we  need  to  consider  how  well  this  approximate  method  performs  in  such  a  setting.  As  a  starting  point  we 
first  consider  simpler  single  particle  observables.  One  observable  with  immediate  relevance  to  AMO 
experiments  is  the  ‘contrast’  or  amplitude  of  the  oscillations  in  a  Ramsey  experiment,  which  for  our  initial  state  is 
given  by  Sx  =  ^  of .  In  figure  1  we  compare  the  time -evolution  of  (Sx )  to  exact  solutions.  For  Ising 

interactions,  Jtj~  =  0  in  equation  ( 1 ),  exact  analytical  expressions  for  the  dynamics  exists  [30, 3 1, 64] .  We  can 
thus  compare  our  DTWA  results  in  a  large  31x31  system.  In  contrast,  for  XT  interactions  (7|  =  0  in 
equation  ( 1)),  no  analytical  solution  is  known  in  2D.  Therefore,  in  this  case  we  resort  to  comparisons  with  ED 
calculation  which  are  limited  to  small  systems  sizes.  In  this  case,  we  choose  a  4  x  5  lattice. 

To  cover  different  regimes,  in  figure  1  we  consider  the  Ising  (panels  (a),  (c))  and  the  XT  (panels  (b),  (d)) 
cases  with  two  different  power-law  decay  exponent,  a  =  3  (panels  (a),  (b)),  and  a  —  l  (panels  (c),  (d)). 
Remarkably,  the  Ising  dynamics  is  exactly  covered  by  the  DTWA  approximation,  an  agreement  that  can  be 
rigorously  justified  [21].  For  the  XT  case  we  also  find  excellent  agreement.  While  for  a  —  3  small  numerical 
differences  are  visible,  for  a  —  1,  the  different  curves  are  nearly  indistinguishable. 
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Figure  1.  Time -evolution  of  the  total  spin  in  the  ^-direction  (or  Ramsey  contrast),  Sx  =  (Sx ).  Compared  are  exact  solutions  (points) 
and  DTWA  results  (lines),  (a),  (c)  Ising  interactions  on  a  3 1  x  31  lattice  with  a  =  3  (panel  (a))  and  a  =  1  (panel  (c)).  (b),  (d)  XY 
interactions  on  a4  x  5  lattice  with  a  =  3  (panel  (b))  and  a  =  1  (panel  (d)). 


4.2.  Spatial  correlations 

We  now  turn  to  checking  the  capability  of  the  DTWA  to  describe  the  time  evolution  of  spatial  two-point 
correlations.  We  again  first  consider  the  case  of  Ising  interactions  where  we  can  compare  the  DTWA  to  an  exact 
analytical  solution  in  large  systems  [3 1 , 32] .  In  particular  we  consider  a  3 1  x  3 1  square  lattice  geometry.  We  study 
the  time- dependence  of  connected  correlation  functions  between  sites  n  and  m; 

=  (GnGm)  ~  (°n  )  (°m)  P  ^  {*»  7>  4,  (9) 

we  calculate  the  correlations  from  the  central  spin  of  the  system,  rn  =  (16,  16),  along  they- direction, 
rm  =  (16,  16  +  j)  (note  that  results  along  the  x- direction  are  identical). 

The  comparisons  for  Ising  interactions  are  summarized  in  figure  2.  We  find  that  in  particular  for 
substantially  long- range  interactions  (case  a  =  1)  the  DTWA  gives  excellent  results.  Only  for  j  ~  1,2  there  are 
small  quantitative  differences  (see  below).  For  shorter  range  interactions  (case  a  =  3)  the  correlations  between 
long-distance  spins  remain  essentially  zero.  The  relevant  short-distance  correlations  are  reproduced  by  the 
DTWA  with  small  deviations  that  are  comparable  to  those  ones  seen  in  the  long-range  interacting  cases.  Overall, 
even  for  a  =  3,  the  spreading  of  the  correlations  is  very  well  (qualitatively)  reproduced  by  the  DTWA.  Note  that 
the  fluctuations  around  zero  for  the  DTWA  solution  are  statistical  because  of  the  finite  number  of  trajectories, 
which  is  of  order  O  ( 105)  in  all  our  calculations. 

By  comparing  the  exact  analytical  solution  of  the  problem  with  the  DTWA  prediction  one  finds  [21]  that 
(of1  of)  (t)  dtwa  =  (ofof)  (t)exact  cos2  (2  tjf ) .  For  our  particular  problem  (note  that  (of )  =  0  for  all  times),  this 
implies  that  the  relative  error  in  the  Cff  correlations  can  be  quantified  as 

1  —  cos2  (  2t/^  )  | .  (10) 

Since  J|  decays  as  a  power-law  with  the  distance  between  spins  i  and  j,  for  short  times,  t  <  7-1,  the  relative  error 
decreases  with  distance  as  well.  In  particular  for  tj  1  it  follows  that  efy  oc  (tj)2j~2a.  Note  that  equation  (10) 
also  implies  that  efy  only  depends  on  the  coupling- strength  between  the  two  spins  in  consideration. 

For  the  XY  model  we  first  compare  the  DTWA  against  ED  results  in  a  small  4x5  lattice.  Due  to  the  small 
system  size,  in  order  to  observe  any  amount  of  linear  spreading  of  correlations  we  have  to  calculate  the 
correlations  from  the  corner  of  the  system,  which  leads  to  additional  boundary  effects.  Specifically,  in  figure  3  we 
calculate  the  correlations  between  the  site  with  rio=  ( 1 ,  1 )  and  sites;  along  the y- direction  with  coordinates 
Tj  =  ( 1 ,  1  +  j).  In  this  XY  case  we  also  find  that  the  DTWA  works  (except  for  deviations  at  the  edge) 
impressively  well,  in  particular  for  C\ J  correlations.  Most  importantly  we  find  that  although  some  oscillation 
seem  not  to  be  well  reproduced,  the  DTWA  accurately  captures  the  spreading  of  the  correlations  and  the  shape 


?yy 


Qyy-> exact  _  QyyX 

lC->Jy  lC>]y 

Qyy,  exact 
Wy 


5 


IOP  Publishing 


NewJ.  Phys.  17  (2015)  065009 


J  Schachenmayer  etal 


Figure  2.  Ising  interaction  benchmark:  time-evolution  of  connected  correlation  functions  Cxxj^  in  a  3 1  x  3 1  lattice  between  the  center 
site  Yic  =  (16,  16)  and  the  sites  rj  =  (16,  16  +  j).  Compared  are  exact  solutions  (upper  plots  in  each  panel)  and  DTWA  results 
(lower  plots),  (a)  a  =  1,  (b)  a  =  3.  The  left  and  right  column  depicts  and  Cj^  ,  respectively. 


of  the  ‘light-cone5  boundary  as  seen  in  figure  3(b).  The  ‘light-cone5  boundary  at  a  threshold  value  Cthres  is 
visualized  in  figure  3  by  a  contour  plot.  Physically  it  corresponds  to  the  propagation  time  required  to  reach  a 
correlation  value  of  Cthres  between  two  spins  separated  by  a  distance;.  Here  we  set  Cthres  —  0-05. 

In  contrast  to  the  Ising  case  where  the  propagation  of  correlations  barely  follows  a  light-cone  spread,  in  the 
XY  model  the  situation  is  more  interesting.  While  for  a  =  1,  the  correlations  build  up  throughout  the  system 
almost  instantaneously,  there  is  a  clear  change  in  behavior  when  going  to  shorter  range  interactions.  In  the  case 
of  a  =  3  a  light- cone  is  expected  to  emerge  [39] ,  and  is  even  visible  in  the  small  system  calculation  shown  in 
figure  3(b).  However,  in  such  a  small  system  boundary  effects  are  expected  to  play  an  important  role.  We  point 
out  that  the  speed  of  the  propagation  of  correlations  in  figure  3  is  essentially  identical  for  C Y-  and  C\ J.  However, 
the  DTWA  result  for  C-j ,  shows  slightly  larger  discrepancies  from  the  exact  solution  than  C{J.  We  also  checked 
the  evolution  of  Cfj  (in  the  XT  case),  which  exhibits  the  same  type  of  correlation  spreading  and  is  equally  well 
reproduced  by  the  DTWA.  However,  in  our  case  this  particular  correlation  is  much  smaller  in  magnitude  than 
and  features  additional  oscillations,  which  is  the  motivation  to  use  in  the  remainder  of  this  article. 

Since  we  will  be  interested  in  light- cones  defined  by  certain  thresholds  of  correlation  values,  one  has  to 
carefully  consider  the  spatial  distribution  of  the  errors  as  they  could  lead  to  wrong  conclusions.  To  test  this 
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Figure  3.  XY  interaction  benchmark:  time -evolution  of  connected  correlation  functions  in  a  4  x  5  lattice  between  the  edge  site 

r ’j0  =  ( 1 ,  1 )  and  the  sites  =(1,1+  j).  Compared  are  exact  diagonalization  solutions  (upper  plots  in  each  panel)  and  DTWA 
results  (lower  plots),  (a)  a  =  1,  (b)  a  =  3.  The  left  and  right  column  is  for  and  CyX  ,  respectively.  White  lines  indicate  contours 

where  C^fy  exceeds  a  threshold  value  of  Cthres  =  0.05. 


dependence  we  compare  the  results  in  a  ID  system  with  1x31  lattice  sites,  where  the  correlations  are  calculated 
from  the  center  [  r y  =  ( 1 ,  1 6)  ] .  In  this  case,  exact  correlations  can  be  easily  computed  by  means  of  t-DMRG 
techniques  [22-25] .  To  average  out  statistical  noise,  we  fit  a  contour  to  a  power  law  of  the  form  Tcthres  cx  f.  The 
error  is  then  calculated  as  the  difference  between  the  DTWA  and  the  exact  t-DMRG  solution  fits  via 
AxJy  =  Results  for  various  threshold  values  and  ranges  of  interactions  are  shown  in  figure  4.  We 

find  that  there  are  small  errors  with  different  sign  at  different  ranges;.  In  case  of  short-range  interactions,  for 
short  distance  correlations  (small;),  the  DTWA  predicts  a  slightly  faster  growth  of  correlations,  while  for  long 
distance  correlations  (large;)  it  tends  to  predict  a  slightly  slower  growth.  For  very  long- ranged  interactions  the 
error  becomes  very  small  and  homogeneously  distributed.  Note  that  in  the  limit  of  nearly  all-to-all  interactions, 
a  1,  the  XY  and  the  Ising  models  become  equivalent  for  fully  symmetric  initial  states.  This  is  because  the  XY 

A  2  A  2 

Hamiltonian  becomes  a  collective  spin- Hamiltonian  ocSx  +  Sy ,  whose  dynamics  is  the  same  as  that  of  the  Ising 

-»2 

model  due  to  the  conservation  of  the  total  collective  spin  S  .  In  this  regime  the  correlations  are  spatially 
homogeneous  and  the  error  is  proportional  to  ejy  oc  ( tj )2. 

In  conclusion,  the  behavior  of  AtJy  with  distance  will  lead  to  small  errors  in  the  predicted  power  law- 
exponent  rj ,  which  become  larger  when  the  range  of  the  interactions  becomes  shorter.  This  is  also  what  we 
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Figure  4.  Error  analysis.  Spatial  dependence  of  the  time  difference  in  fitted  contour  lines.  Different  colors  correspond  to  different 
threshold  values  Cthree  =  0.01,0.02,0.03,0.04,0.05  (from  light  to  dark  color).  We  compare  results  in  a  ID  1  x  3 1  lattice  for  which  we 
can  solve  the  dynamics  exactly  via  t-DMRG  methods.  Correlations  are  calculated  between  sites  with  r \c  =  ( 1 ,  1 6)  and  the  sites 
rj  =  (1,  16  +  j).  In  panels  (a)-(d)  the  interaction  range  increases,  a  =  3,2, 1, 0.5,  respectively. 

Figure  5.  Spreading  of  correlations  in  a  large  3 1  x  31  system  with  XY interactions,  (a)  a  =  4,  (b)  a  =  3,  (c)  a  =  2,and(d) 
a  =  1.  Points  where  the  correlations  exceed  certain  threshold  correlations  Cthres  =  0.01 ,  0.025,  0.05  indicate  a  light-cone  and  are 
shown  as  contour  lines. 


observe  in  figure  7.  In  general,  there  we  find  that  the  various  values  of  rj  are  still  excellently  reproduced  for  a  wide 
range  of  interaction  decay  exponents,  a.  This  confirms  the  validity  of  the  DTWA  to  capture  the  propagation  of 
correlations  in  the  XY  model. 

5.  DTWA  predictions  for  spreading  of  correlations  in  large  systems 

After  having  validated  the  DTWA  method,  we  now  use  this  technique  to  calculate  dynamics  in  a  large  2D  XY 
model.  As  in  our  previous  examples  for  the  Ising  case,  we  focus  on  a  3 1  x  3 1  square  lattice  geometry  and  study 
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Figure  6.  Light  cone  and  contrast  dynamics  for  XY interactions.  Left  panels:  light-cones  indicated  by  contour  lines  for  a  threshold 
correlation  of  Cthres  =  0-05.  Right  panels:  decay  of  Sx  =  (Sx).  (a)  A  large  2D  3 1  x  3 1  system  features  a  rapid  transition  from  light-cone 
behavior  to  all-to  all  physics  at  a  ~  2.  Also  a  qualitative  change  of  behavior  of  the  Sx  decay  occurs  at  a  ~  2.  Left:  power  law  fits, 

Tcthres  a  P  are  shown  as  dotted  lines.  Right:  dashed  lines  show  an  analytical  approximation  to  Sx  (see  text),  (b)  Similar  physics  emerges 
inalDlx31  system  with  a  change  in  behavior  of  the  Sx  decay  at  a  ~  1.  Both  exact  t-DMRG  (points)  and  DTWA  (lines)  are  shown 
and  consistent. 


the  evolution  of  Cy^  .  In  figure  5  we  show  large  system  results  for  decay  exponents  a  =  1,2,3,  and  4.  Again,  by 

defining  the  time  where  the  correlations  exceed  a  certain  value  Cthres> we  can  define  contours  that  indicate  light- 
cones  (see  figure  5  with  Cthres  —  0.01 ,  0.025,  and  0.05).  In  general  we  observe  a  drastic  change  in  the  dynamics  of 
correlations  as  a  is  varied.  While  in  the  case  of  short-range  interactions,  a  =  3,  4,  we  see  a  clear  light- cone-like 
propagation  behavior,  for  a  =  1 ,  2,  this  behavior  breaks  down  rapidly  and  instead  almost  instantaneous 
propagation  of  correlations  is  observed.  This  is  summarized  in  figure  6(a)  (left  panel)  where  we  show  the  light- 
cone  boundary  for  Cthres  =  0-05.  In  the  case  of  a  =  3,  4,  we  can  for  example  easily  fit  a  power-law  curve  to  the 
contour  of  the  form  Tcthree  oc  jn,  with  a  fixed  exponent  rj. 

Interestingly,  the  rapid  change  in  propagation  behavior  is  also  directly  reflected  in  the  time -evolution  of  the 
experimentally  much  more  accessible  observable  SX)  as  demonstrated  in  figure  6(a).  Linked  to  the  disappearance 
of  the  light-cone  (figure  6  left  panel)  at  a  =  2,  the  behavior  of  the  contrast  decay  as  a  function  of  time 
(figure  6(a),  right  panel)  changes.  For  a  =  3,4  after  an  initial  quadratic  decay,  ( Sx )  decays  slowly  (remarkably 
more  slowly  for  a  =  3  than  for  a  =  4).  For  a  ^  2,  however,  these  two  time- scales  disappear  and  (Sx)  exhibits  a 
qualitatively  different  decay.  In  figure  6(b),  we  show  results  for  the  same  calculation  inalDlx31  lattice,  and  a 
corresponding  comparison  to  exact  t-DMRG  calculations.  We  observe  the  same  qualitative  change  in  behavior 
in  the  decay  of  (Sx ),  now  when  crossing  a  ~  1. 

In  the  2D  case,  the  behavior  for  a  ^  2  can  be  understood  semi- analytically.  For  such  long-range 
interactions,  as  explained  above,  we  can  approximately  replace  the  couplings  by  a  constant:  L  «  /eff ,  and  map 

A  2  .  / 

the  dynamics  to  an  Ising  Hamiltonian  —  /eff  Sz  .  The  effective  coupling  constant  /eff  in  our  finite  square  lattice  (M 
spins  in  total)  can  be  determined,  for  example,  by  requiring  that  the  total  energy  of  the  central  spin  interacting 
with  all  other  spins  with  couplings  Jy  is  the  same  as  with  coupling  /eff .  From  the  solution  of  the  Ising  model,  one 
thus  obtains  Sx  ( t)/M  =  cosM_  1  (2t/eff )  which  is  shown  as  dashed  lines  in  figure  6(a)  on  the  right.  Given  that 
there  is  some  arbitrariness  in  defining  a  precise  value  for  /eff ,  we  see  that  the  contrast  decay  for  long-range 
interactions  is  fairly  captured  by  this  simple  model. 

6.  Correlation  dynamics  crossover  with  increasing  range  of  interactions 

Although  the  physics  in  a  ID  chain  seems  to  be  relatively  similar  to  the  2D  system,  a  careful  examination  of 
figure  6  reveals  important  differences.  While  in  the  2D  case,  for  a  <  2  the  contours  are  almost  flat,  in  the  ID  case 
the  contours  still  seems  to  exhibit  a  finite  rj  (recall  Tcthres  oc  jn )  even  for  a  <  1  (i.e.  for  interactions  decaying 
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Figure  7.  Correlation  dynamics  crossover,  (a)  Plots  of  the  contour  lines  reveal  a  power  law  behavior,  Tcthres  oc  t*7  (as  seen  on  a  log-log 
scale).  Points:  DTWA,  solid  line:  power-law  fit,  dashed  line:  t-DMRG  result.  In  2D  a  clear  power  law  is  visible  for  distances  j  >  2, 
which  we  use  as  range  for  the  fit  (in  1 D  we  use  the  whole  range) .  The  power  law  exponent,  r\ ,  as  function  of  the  range  of  interactions,  a, 
is  shown  in  panel  (b)  for  ID  and  (c)  for  2D. 


slower  than  the  system’s  dimensionality).  To  quantify  this  observation  we  systematically  evaluate  rj  as  a  function 
of  the  range  of  inter  actions.  Explicitly  we  vary  a  between  0.5  <  a  <  3  in  ID  and  1.5  <  a  <  4  in  2D  and  set  the 

threshold  value  to  be  Cthres  —  0-05. 

In  figure  7(a)  we  show  selected  examples  of  the  light- cone  contours  for  short  and  long-range  interactions  in 
ID  and  2D  on  a  log-log  scale  (note  that  in  the  ID  case  we  again  compare  our  estimations  to  exact  t-DMRG 
results).  An  interesting  feature  that  we  observe  is  that  in  the  2D  case  the  contour  exhibits  a  clear  power  law 
behavior  only  for  separation  j  >  2.  We  exclude  the  short  distance  correlations  j  <  2  to  perform  the  linear  fit  on 
a  log-log  scale  to  avoid  this  issue. 

In  ID  (see  figure  7  (b))  the  DTWA  nicely  reproduces  the  same  dependance  of  rj  versus  a  seen  in  the  exact 
t-DMRG  calculations  (green  dashed  line).  Although  it  slightly  quantitatively  over-estimates  the  light- cone 
exponent,  it  shows  the  correct  smooth  increase  of  rj  from  zero  to  rj  ~  1  with  increasing  a.  The  situation  in  2D  is 
strikingly  different  and  instead  a  rich  complicated  behavior  is  observed.  Although  statistical  noise  leads  to  non- 
negligible  error-bars,  three  clear  conclusion  can  be  drawn:  (i)  in  contrast  to  ID,  for  a  <  2  (i.e.  for  interactions 
decaying  slower  than  the  system’s  dimensionality)  the  power-law  exponent  of  the  light  cone  is  consistent  with 
rj  =  0;  (ii)  there  is  a  sharp  increase  of  rj  at  a  =  3;  and  (iii)  at  a  ~  4  the  light-cone  behavior  is  consistent  with  a 
linear  causal  region  ( rj  ~  1). 

Given  the  good  agreement  with  exact  calculations  in  the  ID  case,  we  believe  that  the  DTWA  predictions  are 
reliable.  For  the  scenario  in  consideration  (2D  XY  model  with  large  number  of  spins),  no  exact  analytical  or 
numerical  solution  is  available,  and  ultimately  experiments  need  to  provide  a  definite  answer. 

7.  Conclusion  and  outlook 

We  have  used  a  new  numerical  technique,  the  DTWA,  to  study  the  propagation  of  correlations  in  large  2D  XY 
spin  models  with  long-range  interactions,  in  regimes  accessible  to  current  state-of-the  art  experiments  with 


10 


IOP  Publishing 


NewJ.  Phys.  17  (2015)  065009 


J  Schachenmayer  etal 


polar  molecules  or  trapped  ions.  We  benchmarked  this  new  method  in  exactly  solvable  limits  (Ising  interactions 
and  small  systems)  and  found  excellent  agreement.  In  large  systems,  our  method  predicts  a  sharp  change  in  the 
dynamics  exhibited  by  two-point  correlation  and  the  Ramsey  contrast  when  the  decay  exponent  of  the 
interactions  a  crosses  a  —  2.  While  for  a  >  2  a  power- law  light- cone  appears,  the  DTWA  shows  an  additional 
jump  in  the  propagation  speed  of  correlations  at  a  =  3.  For  interactions  with  a  ~  4  the  DTWA  predicts  almost 
linear  light-cone  behavior.  In  the  ID  case  a  power  law  light-cone  is  already  seen  at  decay  exponents  as  low  as 
a  =  0.5,  and  a  nearly  linear  behavior  as  a  ~  3.  We  gained  confidence  in  our  DTWA  prediction  by  direct 
comparisons  to  exact  t-DMRG  calculations  in  ID. 

In  the  future  it  will  be  interesting  to  study  the  nature  of  this  sharp  transition,  not  only  in  2D,  but  also  in  3D. 
This  is  a  regime  currently  accessible  with  polar  molecule  experiments  that  encode  the  spin  degree  of  freedom  in 
rotational  states  coupled  by  dipolar  interactions.  In  such  setups  sharp  changes  in  the  speed  of  correlation 
propagations  could  be  observable.  In  this  implementation  it  will  be  intriguing  to  investigate  the  role  played  by 
the  anisotropic  character  of  the  interactions  and  the  finite  filling  fraction  on  the  light-cone  dynamics.  Systems 
where  retardation  effects  in  the  dipolar  interaction  become  relevant  (e.g.  with  atoms  in  two  electronic  states 
[20, 65] )  could  also  become  excellent  laboratories  for  the  observation  of  DTWA  predictions.  In  many 
implementations  of  spin  models,  dissipation  effects  (due  to  for  example  spontaneous  emission  or  cooperative 
radiation)  compete  with  the  pure  Hamiltonian  dynamics.  In  order  to  model  these  experimentally  relevant 
situations  it  will  be  important  to  adapt  our  technique  to  a  master  equation  formulation  instead  of  pure 
Hamiltonian  evolution. 
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